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Abstract. Toeplitz matrices occur in many mathematical as well as scientific and 
engineering investigations. This paper considers the spectra of banded Toeplitz and quasi- 
Toeplitz matrices with emphasis on non-normal matrices of arbitrarily large order and rel- 
atively small bandwidth. These are the type of matrices that appear in the investigation 
of stability and convergence of difference approximations to partial differential equations. 
Quasi-Toeplitz matrices are the result of non-Dirichlet boundary conditions for the dif- 
ference approximations. The eigenvalue problem for a banded Toeplitz or quasi-Toeplitz 
matrix of large order is, in general, analytically intractable and (for non-normal matri- 
ces) numerically unreliable. An asymptotic (matrix order approaches infinity) approach 
partitions the eigenvalue analysis of a quasi-Toeplitz matrix into two parts, namely the 
analysis for the boundary condition independent spectrum and the analysis for the bound- 
ary condition dependent spectrum. The boundary condition independent spectrum is the 
same as the pure Toeplitz matrix spectrum. Algorithms for computing both parts of the 
spectrum are presented. Examples are used to demonstrate the utility of the algorithms, 
to present some interesting spectra, and to point out some of the numerical difficulties 
encountered when conventional matrix eigenvalue routines are employed for non-normal 
matrices of large order. The analysis for the Toeplitz spectrum also leads to a diagonal 
similarity transformation that improves conventional numerical eigenvalue computations. 
Finally, the algorithm for the asymptotic spectrum is extended to the Toeplitz general- 
ized eigenvalue problem which occurs, for example, in the stability analysis of Pade type 
difference approximations to differential equations. 
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1. Introduction. A Toeplitz matrix has the property that the entries are constant 
along diagonals parallel to the main diagonal. If we define the sequence 

(1.1) a — p, oi — p+i, • • • , a 0 , * * " , Qg — i, dq , where a — p , dq ^ 0 

and p and q axe specified positive integers, then the elements of a banded square Toeplitz 
matrix of order J and bandwidth p + q + 1(< J) are given by 

(1.2) dij — dj — j 

if dj-i is a member of the sequence (1.1) and zero otherwise. Since the eigenvalue analysis 
for triangular matrices is trivial, we have restricted our attention to non-triangular matrices 

1 The main results of this paper were presented at the 14th Biennial Conference on Numerical Analysis, 
Dundee, Scotland, June 25 to 28, 1991. 
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(p, q > 0). For example, if we choose p = 1 and q = 2 and denote the banded Toeplitz 
matrix by A then 


(1.3) 


d$ d\ <32 

a_i <2 q CL \ 0,2 

a_i <3 q Gi Cl2 O 


A = 


0 <2_1 CtQ a l a 2 

< 3—1 CLq d\ 

a— i a,Q _ 


where the bandwidth is 4 and the matrix order J is arbitrary. Our interest is in the spectra, 
i.e., the eigenvalues, for banded Toeplitz matrices of arbitrarily large order. 

A banded quasi - Toeplitz matrix is defined to be a banded Toeplitz matrix where there 
are a limited number of row changes constrained as follows: There are at most p altered 
rows among the first p rows and at most q altered rows among the last q rows. Since p 
and q are fixed and J is assumed to be large, there are only a relatively few altered rows. 
Here we use quasi in the conventional sense meaning almost or nearly. For example, a 
quasi-Toeplitz cousin of the Toeplitz matrix (1.3) is 


(1.4) 


bn 

b 12 

b 13 

b 14 

d-\ 

a o 

ax 

a 2 


a-x 

a o 

ax 


A = 


0 

a-x 

<Z0 a l 



C 2 4 

C 2 3 c 22 

C21 


C14 

Cl3 C 12 

on- 


(The advantage of the unusual indexing for the c entries of A will become apparent in 
Section 4.2). We assume that the modified elements axe constants and therefore indepen- 
dent of J. The bandwidth of a quasi-Toeplitz matrix may be larger than the bandwidth 
of its pure Toeplitz cousin, e.g., the quasi-Toeplitz matrix defined by (1.4) has a larger 
bandwidth than its Toeplitz cousin defined by (1.3). 

Another important relative of a banded Toeplitz matrix is its circulant cousin. Each 
row of a circulant matrix is constructed by cycling the previous row forward one element. 
For example, if the first row of a matrix of order J is 

(1.5) [ao> 0) " ' , 0, 0— p, ■ ,ct_2,fl— i]. 

this process defines the circulant cousin of the Toeplitz matrix defined by (1.2). In partic- 


2 




We will also refer to the matrix defined by (1.5) as the circulant cousin of the quasi-Toeplitz 
matrix, e.g., (1.6) is the circulant cousin of (1.4). A circulant matrix is also Toeplitz 
because the entries are constant along diagonals. The eigenvalues (and eigenvectors) of 
any circulant matrix can be determined analytically [1]. The spectrum for a circulant 
banded Toeplitz matrix is a subset of the spectrum for the doubly infinite (order) banded 
Toeplitz matrix [9]. 

Toeplitz matrices are important in mathematics as well as scientific and engineering ap- 
plications [5,4]. Our interest [14] arose from an investigation of stability and convergence of 
difference approximations to initial-boundary- value problems for partial differential equa- 
tions. The analysis of these difference approximations leads to a study of the spectra for 
banded quasi-Toeplitz matrices of large order. The element changes from the pure Toeplitz 
matrix are the result of applying non-Dirichlet boundary conditions to the difference equa- 
tions. The bandwidth of the matrix is determined by the width of the difference stencil 
and is small compared with the order of the matrix. For example, if the method of lines is 
applied to a scalar hyperbolic partial differential equation and if a four point spatial dif- 
ference approximation is used to approximate the spatial derivative, the stability analysis 
matrix is the quasi-Toeplitz matrix: 


(1.7) 


where the parameters a and (3 are introduced by the numerical boundary conditions [14]. 
If we choose a = 6/5 and /? = —4/5 and calculate the spectrum of (1.7) for a sequence of 
increasing matrix orders, we obtain the spectra shown by the open circles in the complex 
plane in Fig. 1.1. (The abscissa and the ordinate of each eigenvalue, respectively, represent 
the real part and the imaginary part of the eigenvalue.) In the stability analysis of a 
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-10 -10 
(c) J = 24 (d) J = 48 


Figure 1.1. Numerically computed spectra for the quasi- Toeplitz matrix (1-7), with 
a = 6/5 and ft = —4/5 , (o symbol) and the analytically computed spectra for the circulant 
Toeplitz cousin ( + symbol) for matrix order J. 

difference method (and for the purposes of this paper) we axe interested in the spectrum 
for large (approaching infinity) matrix order, e.g., Fig. l.ld. For reference we have plotted 
the spectrum for the circulant cousin of (1.7) in Fig. 1.1. The circulant spectrum is 
indicated by the + symbols. 

In Section 2 the eigenvalue problem for a banded Toeplitz or quasi-Toeplitz matrix is 
formulated as a difference equation with appropriate boundary conditions. In general, the 
resulting boundary value problem is analytically intractable. In Section 3 we assume that 
the order of a quasi-Toeplitz matrix is large and use an asymptotic approach to partition 
the eigenvalues into two independent groups. The partition is between boundary condition 
independent eigenvalues and boundary condition dependent eigenvalues. Algorithms for the 
computation of both parts of the spectrum are described in Section 4. These algorithms 
involve the solution of algebraic equations with degree proportional to the bandwidth of the 
matrix. The algorithm for the boundary condition dependent spectrum closely parallels 
the eigenvalue analysis of Godunov and Ryabenkii [2], Kreiss [7] and Gustafsson, Kreiss, 
and Sundstrom [6]. Plots of the spectra of Toeplitz matrices often exhibit surprising 
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geometrical shapes. Some illustrative examples are shown in Section 5. The algorithms for 
computing the spectra are applied to both Toeplitz and quasi-Toeplitz matrices in Section 
6. If conventional numerical eigenvalue algorithms are used to compute the spectra for non- 
normal matrices of large order, errors in the spectra may be encountered. In Section 7 we 
introduce a simple similarity transformation that improves the accuracy of conventional 
eigenvalue algorithms for a large class of banded Toeplitz and quasi-Toeplitz matrices. 
Finally, the Toeplitz eigenvalue algorithm is extended to the generalized eigenvalue problem 
in Section 8. 

2. The eigenvalue problem: a difference equation representation. Let A 

represent a J x J banded Toeplitz matrix defined by (1.2). The eigenvalue problem for A 
is defined by 

(2.1) A <f> = A <f> 

where <f> is the eigenvector and A is the eigenvalue. The indexing of the eigenvector <f> is 
arbitrary and we choose 

(2.2) (f> T = [<f>i,<f>2,” m Aj-iAj]- 

The eigenvalue problem (2.1) is equivalent to the system of homogeneous difference 
equations 

7 

(2.3) ^ j+m = A (f>j, j 1, 2, • ■ ■ , J 

m=—p 

with p homogeneous Dirichlet boundary conditions at the left boundary and q homogeneous 
boundary conditions at the right boundary : 

(2.4a,b) = 0, in — 0 , 1, • ’ • ,p 1, j-\-m 1) i Q- 

Note that fictitious points have been introduced in defining the boundary conditions (2.4). 
It is apparent that (2.3) is the generic equation of the system (2.1). 

Each of the equations (2.3) is a (p + g)-th order difference equation for <fij and we look 
for a solution of the form 

(2.5) 4>j = /c J 

where k is a complex constant. Insertion of (2.5) into (2.3) yields 

7 

(2.6) A= 

771 — — p 

This is an algebraic equation of degree p + q in the unknown k and in the following analysis 
we assume the roots are distinct. The general solution of (2.3) is 

p+q 

(2.7) <f>i = fim K in 

m= 1 
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where the K m 's axe the roots of (2.6) and the /? m ’s are arbitrary constants. The constants 
/? m are determined by inserting (2.7) into the boundary conditions (2.4). If the bandwidth 
is three, i.e., p = q = 1, the eigenvalue problem for the pure Toeplitz problem can be 
solved analytically [10,15]. However, if the bandwidth is greater than three, the problem 
is analytically intractable. 

For a quasi- Toeplitz matrix the difference equation formulation of the eigenvalue prob- 
lem differs from the pure Toeplitz problem only in the boundary conditions. The difference 
equations (2.3) are still valid if appropriate extrapolation boundary conditions are gener- 
ated. These extrapolation boundary conditions depend on the matrix entries that change 
the matrix from Toeplitz to quasi-Toeplitz. Extrapolation boundary conditions have the 
general form 

(2.8a,b) = 0, vn — 0, 1, • • • , p 1, Q m (E^f> j+m = 0> m >9 

where P- m {E) and Q m {E ) are polynomial operators and E is the shift operator defined 
by Ecf>j = <f>j+ 1 . The subscripts on the operators indicate that, in general, the operators 
are different for each value of m. As an example, the extrapolation boundary conditions 
for the matrix (1.7) are 

(2.9a, b,c) P 0 (E)<f>o = 0, Qi(E)<j>j+\ = 0, Q2{E)4>j+2 = 0 

where, after some algebra, one finds 

(2.10a, b, c) Po(E) = (E — l) 3 (3aF — 1), Qi(E) = 1, Qi{E) = (1 — E~ 1 ) 3 (l + 6/3E 1 ). 

Although the extrapolation boundary conditions are easy to derive, we will find it more 
convenient to implement an eigenvalue algorithm which uses the boundary difference equa- 
tions directly. For example, the first row and the last two rows of the matrix (1.4) have 
entry changes. For the eigenvalue problem (2.1) the difference equations (2.3) are still 
valid at j = 2, 3, • ■ • , J — 2 and modified difference equations apply at j = 1, J — 1, J. The 
general form of difference equations for the left boundary is 

(2.11) B^ = AD^ 

- T 

where B is a rectangular matrix with p rows and r(< J) columns and <f> = [<f>\ , <j >2 , ■ - • , 4> r ] ■ 

The matrix D is a p x r rectangular diagonal matrix with unit entries on the diagonal. 
Similar equations apply at the right boundary. 

The eigenvalue problem for the circulant cousin of (1.2) is (2.3) with periodic boundary 
conditions: 

(2.12) <f>j = <j>j+j. 

It is well-known that the eigenvalue problem for any circulant matrix can be solved ana- 
lytically. By inserting (2.5) into (2.12), one obtains 

(2.13) k j = 1. 
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Hence the k’s axe the J roots of unity: 

(2.14) ki = e' 91 , where 9f=2£n/J, £=1,2, •••,</. 

The circulant spectrum is obtained by substituting (2.14) into (2.6): 

(2.15) A* = Y, a m^ im9 ‘- 

771 = —p 

As an example, the spectrum of the circulant matrix (1.6) with coefficients 

(2.16) [a_ 1 ,a 0 ,a 1 ,a 2 ] = [-1/3, -1/2 , 1, -1/6] 

reduces to the compact form: 

(2.17) A i = —4:wj/Z — i[l + 2wt/Z] sin0/, 
where 

(2.18) u>t = sin 2 (0*/2), 9t—2£Tr/J, £=1,2, ■ ■ ■ , J. 

The eigenvalue locus is the oval shaped dashed line plotted in Fig. 5.1. In equation (2.16) 
and the remainder of this paper we use the underline to indicate the main diagonal element. 

3. Partition of the spectrum. If we assume that the order J of a quasi-Toeplitz 
matrix is large (J — ► oo) we can partition the eigenvalues into two groups. The parti- 
tion is between (a) eigenvalues which are independent of the boundary conditions and (b) 
eigenvalues which are dependent on the boundary conditions. Group (b) may be empty. 

We assume without loss of generality that the p + q roots of the algebraic equation (2.6) 
are ordered as 

(3.1) |«l| < |«2 1 < < |*7>| < K+l| < < |«p+?|- 

It should be noted that each k is a function of A. 

The crucial inequality in (3.1) is between |/c p | and |kj,+i|. As J — * oo there are only two 
choices: either 

(3.2a,b) |k p | = |k p +i| or |/c p | < |/c p +i|. 

This simple observation divides the asymptotic spectrum into two groups. In this section, 
we briefly describe the properties of the two groups. Schmidt and Spitzer [9] proved the 
equivalence of the eigenvalues associated with the /c’s of (3.2a) and the asymptotic ( J — ► oo) 
spectrum of the banded Toeplitz matrix. 

We first consider the case of equality (3.2a): 

(3.3) group (a): KI = |« P +i|- 
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For large J (J — * oo), the k’s of equal modulus (3.3) couple the left and right boundaries 
for the eigenfunction given by (2.7). On the other hand, the equality (3.3) allows one to 
find the asymptotic spectrum which is independent of the boundary conditions. This may 
seem like a paradox, but the algorithmic details in the next section clarify how it is possible 
for the k’s and the spectrum to be independent of the boundary conditions. 

Next we consider the inequality (3.2b): 

(3.4) group (b): \k p \ < |k p+1 |. 

In this case, the left and right boundaries are uncoupled for large J (J — ► oo). Furthermore, 
the corresponding eigenvalues are dependent on the boundary conditions. Eigenvalues in 
group (b) depend intimately on the choice of boundary conditions and for certain choices 
there will be no eigenvalues in this group. One choice of boundary conditions which 
introduces no eigenvalues of group (b) is the Dirichlet boundary conditions (2.4) (see, 
e.g., Section 6.2). Since Dirichlet boundary conditions lead to a pure Toeplitz eigenvalue 
problem, we conclude that the spectrum for a pure Toeplitz matrix is composed only of the 
boundary condition independent eigenvalues. 

4. Algorithms for computing the spectra. In this section we present algorithms for 
computing the asymptotic spectra of Toeplitz and quasi-Toeplitz matrices. We separate the 
algorithms according to the two groups of eigenvalues: (a) boundary condition independent 
(pure Toeplitz) and (b) boundary condition dependent. The application of each algorithm 
involves the solution of algebraic equations whose degree is proportional to the bandwidth 
p + q + 1 of the matrix. 

4.1. Boundary condition independent (Toeplitz) spectrum. The algorithm for 
the pure Toeplitz spectrum is based on the observation that separated the spectrum into 
two groups, i.e., we seek those eigensolutions whose k’s satisfy (3.1) with the equality 
(3.3). First we seek the solutions of the algebraic equation (2.6) with two distinct but 
equal modulus k’s. We denote these k’s by 

(4.1) K a = ke 1 ^ 1 , Kfc = ke~'^ 1 , 0 < < n 

where k is a complex variable. For computational convenience, the interval [0, tt] is divided 
into M + 1 subintervals of size A rp where (M + 1)A^ = n and rft = ^A ip, i = 1, 2, ■ • ■ , M. 
(The choice of the integer M determines the number of computed points on the spectrum, 
i.e., the resolution of the spectrum, but does not affect the accuracy of the computed 
points.) An equation for k (independent of A) is easily obtained. Since n a and k& must 
each give the same value of A when inserted in (2.6), we can eliminate A and obtain a 
polynomial in k with coefficients which are functions of the Toeplitz matrix elements and 
x[>f. In particular, 

? 

(4.2) A(k 0 ) ~ 

m= —p 

and since 

(4.3) A(k 0 ) = A(kj) 
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A («,) = £ 

m=— p 


one obtains a polynomial equation of order p + q for k: 

(4.4) 53 °m sin(m^)K m = 0. 

m=— p 

The algorithm for the boundary condition independent spectrum is: 

Algorithm 4.1 

For each of the M values of ifig: 

(i) Solve equation (4.4) for the p + q roots k. 

(ii) Check the corresponding k’s to determine if they satisfy (3.1) with (3.3), i.e., that k„ 
and Kb from (4.1) are in fact k p and k p +j. This is done as follows. 

For each k obtained in (i) 

(iia) Calculate K a and kj from (4.1) and insert K a (or Kb) in (2.6) to find A. Substitute 
this value for A back into (2.6) and solve the algebraic equation to determine the remaining 
(p + q - 2) k ’ s . 

(iib) If the p + q k’s satisfy the inequality (3.1) with (3.3), i.e., |k„| = |k(,| = |k p | = |k p +i|, 
then A is a point on the asymptotic Toeplitz spectrum. If the k’s do not satisfy the 
inequality test, discard the A. 

Return to (iia) until all k’s from (i) have been tested. 

Replace ipg by rpg+\ and return to step (i). 

It should be emphasized that the algorithm (4.1) requires the solution of many algebraic 
equations but their degree is proportional to the bandwidth of the matrix and not the 
order of the matrix. Note that the boundary conditions do not enter the calculation. An 
implementation of the algorithm in MATLAB is given in the Appendix . The algorithm 
will become more transparent in the examples of Section 6. 

4.2. Boundary condition dependent spectrum. The algorithm for the bound- 
ary condition dependent spectrum closely parallels the eigenvalue analysis of Kreiss [7] 
and Gustafsson, Kreiss, and Sundstrom [6]. There are however some differences since 
their analysis for hyperbolic initial-boundary-value problems makes two assumptions: (I) 
Cauchy stability (an assumption about the eigenvalues of the Toeplitz matrix’s circulant 
cousin) and (II) the quarter-plane eigensolutions are unstable (i.e., the eigenvalues are in 
a particular part of the complex plane). Since we are interested in the entire spectrum 
for any quasi-Toeplitz matrix, we do not impose conditions (I) and (II). The relaxation of 
these two conditions leads to a slightly more complicated algorithm. 

The algorithm for the boundary condition dependent part of the spectrum is straight- 
forward. We seek p (or q) values of k which satisfy (2.6) and satisfy the left p (or right 
q) boundary difference equations. The remaining q (or p) k’s must satisfy (3.1) with 
inequality (3.4). The algorithm for the eigenvalues associated with the left boundary is: 
Algorithm 4.2 

(i) Assume a solution of the form 

( 4 - 5 ) = 53 

m=l 
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Substitute <t>j from (4.5) into the p left boundary difference equations (2.11) to obtain p 
equations for the 2p + 1 unknowns: p K m ’s, p /? m ’s, and A. These p equations are linear 
and homogeneous in the 0's therefore the nontrivial solutions reduce the system to a single 
equation for the p k's and A. An additional p equations are required and they are derived 
from the algebraic equation (2.6) since each of the p k's must give the same value of A. 
This system of polynomial equations can be reduced to a single polynomial equation for a 
single variable and we choose k p . 

(ii) The roots of the algebraic equation from (i) contain all candidate k p 's for the (left) 
boundary condition dependent spectrum. For each k p the corresponding K m 's must be 
checked. This is accomplished by calculating the corresponding A (k = k p ) from (2.6) and 
then computing the p + q roots of (2.6) for the given A. Finally, we compare those roots 
with the p K m 's which satisfy the boundary difference equations (2.11). If (3.1) and (3.4) 
are both satisfied, the A is an eigenvalue in the boundary condition dependent part of the 
spectrum. 

The algorithm for the eigenvalues associated with the right boundary is similar. If the 
algorithm has been implemented for the left boundary problem, it can be used for the right 
boundary problem by pre- and post-multiplying A by the permutation matrix 


(4.6) 


0 


P = 


1 


1 1 




1 


1 


O 


i.e., PAP. Note that P -1 = P. The matrix PAP is most easily obtained by interchanging 
the columns of A “left to right” and then interchanging the rows “up and down”. The 
boundary difference equations associated with the left boundary for PAP are 


(4.7) 


= AD^. 


As an example, the matrices C and D associated with PAP where A is (1.4) are 


(4.8) 


Cil Ci2 Cj3 Ci4 

D — 

'1 000' 

C 21 C 2 2 C 23 C 24 _ 


0 1 00 


This explains the somewhat unconventional indexing for the c elements of (1.4). (Note 
that a pure Toeplitz matrix A is persymmetric, i.e., it is symmetric about its northeast- 
southwest diagonal [3]. Consequently A T = PAP. However if A is quasi- Toeplitz then 
PAP is not the transpose of A.) 

The details of the algorithm are illustrated by the examples of Section 6. Note that there 
may be no eigenvalues from the boundary condition dependent part of the spectrum even 
if non-Dirichlet boundary conditions are specified, i.e., the set of eigenvalues corresponding 
to group (b) defined by (3.4) may be empty. 
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The algorithm for the boundary condition dependent spectrum is straightforward but 
finding a polynomial equation in a single variable can be algebraically involved and the 
degree of the polynomial can be large. The system of equations for the k's in (i) are 
symmetric polynomials therefore the algebraic complexity can be somewhat reduced by 
the use of the elementary symmetric functions. We have implemented algorithm 4.2 in 
MACSYMA and it performs satisfactorily on an IRIS workstation for the matrix B with 
two rows and a quasi- Toeplitz matrix A with bandwidths up to five. 

4.3. Alternative boundary condition dependent algorithm. The algebraic com- 
plexity of the boundary condition dependent algorithm of Section 4.2 led us to consider 
alternative methods. In the GKS (Gustafsson, Kreiss and Sundstrom [6]) theory a GKS 
(i.e., unstable ) eigenvalue is a member of the boundary condition dependent eigenvalue set 
of this paper. Kreiss [7] recognized that it is often very difficult to find GKS eigenvalues by 
analytical methods. He proposed finding GKS eigenvalues (if they exist) by a numerical 
eigenvalue computation of the finite-domain problem where the boundary condition (or 
conditions) to be checked are on the left boundary and Dirichlet boundary conditions are 
imposed on the right boundary. Kreiss showed that GKS eigenvalues converge exponen- 
tially fast with increasing J . We found however difficulties in numerically delineating the 
GKS eigenvalues (or other boundary condition dependent eigenvalues) for non- normal ma- 
trices. The numerical difficulties are often alleviated if the matrix is rescaled as described 
in Section 7. 

We have found the following algorithm to be useful although not always completely 
reliable since it will sometimes “miss” an eigenvalue. The reliability is improved if the 
algorithm is repeated for several matrix orders. 

Algorithm 4.3 

(i) Place the boundary conditions of interest on the left boundary and Dirichlet conditions 
on the right boundary. 

(ii) Apply the similarity transformation described in Section 7 and use a conventional 
eigenvalue algorithm to find the eigenvalues for large J. 

(iii) Test all eigenvalues found in (ii) using Part (ii) of algorithm 4.2. 

5. Numerical application of algorithms. Plots of the spectra of Toeplitz matrices 
present some interesting and sometimes surprising geometrical shapes. In this section we 
present some spectra which were obtained with the algorithms described in section 4. The 
matrix coefficients were selected more or less at random to give a flavor of the geometric 
patterns that may occur. 

In each figure we show the spectra plotted in the complex plane. Open circles represent 
numerically computed eigenvalues for the matrix order (J) specified in the figure legend. 
The “solid” line (which is formed by closely spaced dots) is the asymptotic spectrum for 
the Toeplitz matrix and was computed using algorithm 4.1. The star symbol (*) is used to 
represent the boundary condition dependent part of the asymptotic spectrum computed 
using algorithm 4.2 or 4.3. The dashed line represents the spectrum for the circulant cousin 
of the pure Toeplitz matrix. 

5.1. Toeplitz spectra. First we summarize the known properties of the asymptotic 
spectra of pure Toeplitz matrices and show some examples which illustrate these properties. 
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Figure 5.1. Toeplitz matrix [-1/3, -1/2, 1, -1/6] spectra “fill in” for increasing matrix 
order J. The symbol o denotes a numerically computed eigenvalue. Asymptotic spectra 
(algorithm f.l) are denoted by “solid” lines and asymptotic circulant spectra (2.17) by 
dashed lines. 

Let A j denote a J x J pure Toeplitz matrix where A(A j) = {Ai, A 2 , • • • , A j) denotes the 
spectrum. Consider a sequence of matrices A j where the dimension J tends to infinity. 
As J increases a subset of the complex plane is “filled in” by the eigenvalues A(Aj) as 
J — ► oo and one obtains the asymptotic spectrum [9]. This is illustrated in Fig. 5.1 
for a quadridiagonal pure Toeplitz matrix defined by the sequence (2.16). Note that the 
spectrum A(A/) is close to the asymptotic spectrum even for small values of J . 

The asymptotic spectrum of a banded pure Toeplitz matrix is compact and locally 
it consists of analytical arcs [9]. Furthermore, it possesses no isolated points [13]. This 
property will be of interest in considering the asymptotic spectra of quasi-Toeplitz matrices. 
Ullman [13] has shown that the asymptotic spectrum is connected. Asymptotic spectra 
illustrating these properties axe plotted in Figs. 5.2 and 5.3. The defining sequence (1.1) 
is shown under each figure. Algorithm 4.1 is not restricted to Toeplitz matrices with 
real coefficients and this is demonstrated in Fig. 5.3d. As the bandwidth of the matrix 
increases, the asymptotic spectrum becomes increasingly complex. 
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-2 0 2 -2 0 2 


(a) [1,0, 0,1] (b) [.5, 0,0, 1.5] 



(c) [.75, .75,0, .25] ( d) [1, -1,0, .25] 

FIGURE 5.2. Asymptotic spectra for pure Toeplitz matrices (algorithm f.l) indicated by 
“solid” lines and the asymptotic spectra for the circulant cousin Toeplitz matrices indicated 
by dashed lines. 

According to Schmidt and Spitzer [9] it is not easy to give a simple geometric description 
of the asymptotic Toeplitz spectrum for non-normal matrices. An exception is the special 
case where the defining sequence is [a_ p ,0, , 0, a g ], i.e., all members of the sequence 

are zero except the two “outriders”. In this case the spectrum is a star-shaped figure. 
Examples of star-shaped spectra for quadridiagonal matrices axe shown in Figs. 5.2a,b 
and for a septadiagonal matrix in Fig. 5.3a. The spectrum for the tridiagonal case is 
known analytically and consists of a straight line segment (connecting the foci of the 
ellipse which is the spectrum of its circulant cousin) [15]. 

It is a general result [9] that the asymptotic Toeplitz spectrum is “enclosed” by the 
spectrum of its circulant cousin. This is illustrated in Figs. 5.2 and 5.3. However, for 
finite J the spectrum A(Aj) of A is not necessarily enclosed by the asymptotic spectrum 
of the circulant cousin of A. This is illustrated in Fig. 5.4. 

5.2. Quasi- Toeplitz spectra. The properties of the asymptotic spectrum of a quasi- 
Toeplitz matrix are illustrated in Fig. 5.5 for the matrix (1.7). The asymptotic spectrum 
consists of the asymptotic pure Toeplitz spectrum, group (a), plus (possible) isolated eigen- 
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1 
5 
0 
5 
1 

-0.5 0 0.5 1 

(a) [1,0, 0,0, 0,0, .75] (6) [.45, -.45, 0, -.05, -.05] 

4 
2 
0 
-2 
-4 

-10 12 -10 1 
(c) [—.2, —.2, —.2,^1, —.2, .2, —.4] (<f) [.7i, .3, 0, .5i, .3 + i] 

FIGURE 5.3. Asymptotic spectra for Toeplitz matrices (algorithm 4-1) indicated by “solid” 
lines and their circulant cousin Toeplitz matrices indicated by dashed lines. 

values corresponding to boundary condition dependent eigenvalues, group (b). The isolated 
eigenvalues (see, e.g., Fig. 5.5b) must be boundary condition dependent eigenvalues as- 
sociated with group (b) because the pure Toeplitz spectrum can have no isolated points 
[9]. Note also that, in contrast to the pure Toeplitz spectrum, the quasi- Toeplitz spectrum 
may have (isolated) eigenvalues outside the circulant spectrum, e.g., Figs. 5.5b,c,d. For 
a = /? = 0 the matrix (1.7) has no boundary condition dependent eigenvalues and the 
asymptotic quasi- Toeplitz spectrum is the pure Toeplitz spectrum shown in Fig. 5.5a. In 
this example group (b) is empty. A fundamental property of the eigenvalues associated 
with group (b) is that the “left” and “right” boundary dependent eigenvalues are uncou- 
pled. This is illustrated in Fig. 5.5b,c,d. For (3 = — .8 there axe three isolated eigenvalues 
to the right of the pure Toeplitz spectrum. These eigenvalues axe completely independent 
of the parameter a of the left boundary condition (see matrix (1.7)). This fact is illustrated 
by comparing Figs. 5.5b,c where /? is fixed and a is changed. For a = 1.2 there is one 
isolated eigenvalue to the left of the circulant spectrum. This eigenvalue is independent of 
the parameter (3 of the right boundary condition as is indicated by comparing Figs. 5.5b,d. 

6. Analytical application of algorithms. If the bandwidth of the matrix is suffi- 
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(c) J = 24 (d) J = 48 

FIGURE 5.4. Spectra for finite matrix order J (o symbol) and the asymptotic circu- 
lant spectrum (dashed lines). Asymptotic Toeplitz spectra indicated by “solid” lines. The 
Toeplitz matrix is [.7, .1,0, .4, .6]. 

ciently small (sa 3) the algorithms of Section 4 can be used as analytical methods. In this 
section we elucidate the algorithms of the previous section by considering small bandwidth 
examples for both Toeplitz and quasi-Toeplitz matrices. 

6.1. Toeplitz tridiagonal matrix. A tridiagonal Toeplitz matrix is defined by the 
sequence [a_i, ao, ai] where p — q = 1 (see (1.1)). For this case the spectrum for arbitrary 
order J can be determined analytically. The solution is “well-known” although the eigen- 
value formula does not appear in many references. A derivation is given by Smith [10, p. 
**]. The algebraic equation (2.6) becomes 

(6.1) A = a_i k~ 1 + ao + «i «• 

By rewriting (6.1), one obtains 

(6.2) a\K 2 + (ao — A)/c + a_i = 0. 
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FIGURE 5.5. Asymptotic spectra for the quasi- Toeplitz matrix (1-7). Boundary condition 
independent spectra computed with algorithm f.i (‘‘solid” line) and boundary condition 
dependent spectra (* symbol) computed with algorithm 4-2. 

We denote the roots of this quadratic by K\,k 2 . The product of the roots is 

a_i 

(6.3) kik 2 . 

v ■ «i 

We follow the algorithm for the boundary condition independent part of the spectrum 
from Section 4.1. Equation (4.4) for k reduces to 

(6.4) a\k 2 — a_i = 0. 

In step (i) we solve for the roots of (6.4), i.e. , 

(6.5) k = ±\/a_i/ax. 

In step (iia) we find 

(6.6a,b) K a = \J a-i/aie'^ , Kb = 'a~\Ja[e X ^ 1 
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where we have chosen the positive square root. (The negative square root gives the same 
spectrum.) Since there are only two roots of (6.1), i.e., and « 2 j they must be K a 
and /tfr. They obviously have equal modulus, therefore equality (3.3) is satisfied and the 
corresponding A is part of the Toeplitz spectrum. (In this simple example step (ii) of the 
algorithm 4.1 is rather trivial.) Substitution of K a given by (6.6a) into (6.1) yields 

(6.7) A* = ao + 2-^/aid-i cos i/>/, 0 < ipt < tt 

which is the asymptotic spectrum for a tridiagonal Toeplitz matrix. 

Although not part of the asymptotic analysis we note that for the tridiagonal case, the 
roots Kx and «2 have equal modulus regardless of the order J of the matrix. In fact, the 
exact roots are given by 

(6.8) K\ = ke'^ 1 , «2 = ke ~ 1 ^ 1 , where rpt = £ir/(J + 1), £ = 1,2 ,•••,</ 

and k is given by (6.5) (see, e.g. [15]). Consequently, the exact spectrum for a J x J 
matrix is (6.7) where tpt is defined in (6.8). It should be emphasized that if the bandwidth 
of a Toeplitz matrix is larger than three, then algorithm 4.1 is, in general, exact only in 
the asymptotic limit J — * oo. 

6.2. Quasi- Toeplitz tridiagonal matrix. The asymptotic spectrum for a quasi- 
Toeplitz tridiagonal matrix is given by (6.7) plus any boundary condition dependent eigen- 
values. We first show that a pure Toeplitz tridiagonal matrix does not have any boundary 
condition dependent eigenvalues, i.e., the homogeneous (Dirichlet) boundary conditions 
(2.4) introduce no boundary condition dependent eigenvalues. For a tridiagonal matrix 
the homogeneous conditions (2.4) are simply 

(6.9a, b) <f>o = 0, (f>j + 1 = 0. 

By substituting (4.5), i.e., 

(6.10) <f>j = (3x4 

into (6.9a), one concludes that either (3\ = 0 or «i = 0. Since the product of the roots 
(6.3) is nonzero, then /?i = 0 and there is no nontrivial eigenfunction and no eigenvalue is 
introduced by the (left) boundary condition (6.9a). By substituting (6.10) into (6.9b), one 
finds the similar result that no eigenvalue is introduced by the (right) boundary condition 
(6.9b). 

For the general Toeplitz matrix there are p homogeneous Dirichlet boundary conditions 
at the left boundary and q homogeneous Dirichlet boundary conditions at the right bound- 
ary, i.e., (2.4a,b). Substitution of (4.5) into the left boundary conditions leads to a homo- 
geneous system V/? = 0 where V is a p x p Vandermonde matrix and (3 = [(3\,(3 2 , — • , (3 p \ T . 
The determinant of V equals zero iff two «’s are equal which violates the assumption of 
distinct k's. Therefore Dirichlet boundary conditions introduce no boundary condition 
dependent eigenvalues. 
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Consider now boundary conditions that introduce nontrivial solutions. Suppose that 
when (6.10) is introduced into the left boundary condition we obtain 


(6.11) *1 = *o- 

where K a is some nonzero complex constant. Does this introduce an eigenvalue? It does if 
inequality (3.4) is satisfied, i.e., |«i| < |/c 2 1. From (6.3) one obtains 

(6.12) |k 2 | = |a-i/ai|/|«:«r|. 


Inequality (3.4) is satisfied if 

(6.13) |*i| = 1*0-1 < VV-i/ai| 


and the corresponding eigenvalue A is obtained by substituting into (6.1). 

We give an example of a quasi-Toeplitz matrix where there is a boundary condition 
dependent eigenvalue. Consider the matrix 


(6.14) 


■ 0-2 2 

-10 1 O 


O -10 1 

-1 0 . 


Here 

(6.15) a_ i — —1, ao =0, a x = 1, &n = 0, &12 = —2, &13 = 2. 

The eigenvalue problem is (2.3) with p = q = 1 and the boundary difference equation 
(2.11) is 


(6.16) 


— 2(^2 + 2^3 — A(£i. 


Substitution of (6.10) into (6.16) yields 

(6.17) — 2«i + 2*f = A. 

The second equation relating /cj and A is (6.2) with * = *i and coefficients (6.15): 

(6.18) - A*i — 1 = 0. 

Elimination of A from (6.17) and (6.18) yields 

(6.19) (*1 — 1) 2 (2 *i + 1) = 0. 
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The roots of (6.19) axe 

(6.20) «i = 1, 1, —1/2. 

In the tridiagonal case, step (ii) of algorithm 4.2 is simplified to checking inequality (6.13). 
The inequality (6.13), with (6.15), 


|«i| = |«<r| < 1 

is satisfied for = —1/2. Substitution of ki = —1/2 into (6.17) gives A = 3/2. 

The asymptotic spectrum of the quasi- Toeplitz matrix (6.14) consists of two parts. The 
boundary condition independent (pure Toeplitz) spectrum is (6.7): 

(6.21) A = i2cost/>, 0 < ip < tt 

and the boundary condition dependent part of the spectrum is A = 3/2. Consequently, the 
asymptotic spectrum is a line segment on the imaginary axis between ±2? and the single 
eigenvalue A = 3/2. 

6.3. Toeplitz quadridiagonal matrix. A more interesting example is the quadridi- 
agonal matrix (1.3) where p = 1 and q = 2. The algebraic equation (2.6) is 

(6.22) A = a_i« -1 + ao + aiK + a2/c 2 . 

For the boundary condition independent spectrum, the three roots of the cubic (6.22) must 
satisfy (3.1) and (3.3) 

(6.23) |«i| = [« 2 1 < |k 3 |. 

The product of the roots of (6.22) is — a_i/a 2 and therefore 

(6.24) |«l||«2||K3| = | a -l/ a 2|- 
The polynomial (4.4) for k becomes 

a_ i sin (—xpt)k~ l + aj sin(i/><)/c + 02 sin(2^y)/c 2 = 0 


or 

(6.25) 2a2(cos ipf)k 3 + a\k 2 — a_i = 0. • 

In (i) of algorithm 4.1 we solve (6.25) numerically for specified coefficients at. 

In (iia) of the algorithm 4.1 we proceed as follows. For each root k we substitute 
K a = ke"^ 1 into (6.22) and calculate the corresponding A. Then with A given, we solve the 
cubic (6.22) for the roots /«i,K 2 ,K 3 - Of course, two of the roots are already known, i.e., 
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K a = ke'^ 1 and /cj, = /ce - '^. In (iib) we check to see if the three roots satisfy inequality 
(6.23): 

(6.26) |««| = M = |«i| = M < M- 

If this inequality is satisfied, then the corresponding A is a point of the boundary condition 
independent spectrum. If the k’s fail to satisfy inequality (6.26), we discard the A. We 
repeat (i) and (ii) for all of the t/Vs where l = 1, 2, • • • , M. 


For the particular example under consideration, the test (6.26) can be simplified as 
follows. Recall that the roots of equal modulus K\ and n 2 are denoted by (4.1) and rewrite 
(6.23) as 

(6.27) |k| = |«| < |«3 1 
and (6.24) as 

(6.28) |«|l«||/c 3 | = |a_i/a 2 |- 
Inequality (6.27) is satisfied if 

(6.29) \k\ < |a-i/a 2 [ 1/3 . 

Hence if |k| satisfies (6.29), then the corresponding A computed by inserting K a = ke'^ 
into (6.22) is a point of the boundary condition independent spectrum. 

The spectrum for the Toeplitz matrix (1.3) with coefficients (2.16) was computed with 
algorithm 4.1 (see also the Appendix) and is the “solid” line in each of the figures of Fig. 
5.1. The spectrum for the circulant cousin (1.6) is the dashed line in each of the figures of 
Fig. 5.1 (see (2.17)). 

6.4. Quasi- Toeplitz quadridiagonal matrix. An example of a quasi-Toeplitz 
quadri diagonal matrix is (1.4) where p = 1 and q = 2. The eigenvalue problem is equivalent 
to (2.3), i.e., 

(6.30) \<f>j = a-i4>j-i + a 0 <f)j + ai<£>+i + a 2 <f>j+ 2 i j = 2,3, •••,«/ 
with left boundary difference equation (2.11): 

(6.31) hl<f>l + b\ 2 (f> 2 + i>i3<^>3 + bi4(f>4 = 

The spectrum for the matrix (1.4) is the spectrum for the pure Toeplitz matrix (plotted 
in Fig. 5.5b for coefficients (2.16)) plus any boundary condition dependent eigenvalues. 
For the boundary condition dependent spectrum the roots of the cubic (6.22) must satisfy 

(6.32) |«i| < |« 2 | < |«s| 
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(see (3.1) and (3.4)). 

In the analysis of the boundary condition dependent part of the spectrum, the left and 
right boundary conditions are uncoupled and consequently we consider each boundary 
separately. 

6.4.1. Left boundary dependent spectrum. The left boundary difference equation 
is given by (6.31). We follow the algorithm for the boundary condition dependent part of 
the spectrum given in Section 4.2. Recall that for this example p = 1 and in the first step 
(i) we look for a solution of the form (4.5) 

(6.33) <f>j = 0i k{. 

Insertion of this equation into the boundary difference equation (6.31) yields 

(6.34) b\ i + b\2ft\ + f>i3 K i + &i4 K i = A. 

As a particular example we consider the special case for the matrix (1.4) where the 
coefficients at are given by (2.16) and the coefficients b\j are given by 

(6.35) &ii — — ct — 3/2, b \ 2 = 3a 4- 2, b \ 3 = — 3a — 1/2, 614 — a. 

The algebraic equation (6.22) with k = K\ becomes 

2 

(6.36) 6A — 3 4- 6/ci — k\ = —(/ci — 1)(/Cj — 5/ci — 2)//ci 

and the boundary difference equation (6.34) with coefficients (6.35) becomes 

(6.37) -(a + 3/2) + (3a 4- 2)/c x - (3a 4- l/2)/c? 4- a/c? = A. 

If we eliminate A from (6.36) and (6.37), the polynomial for k\ is 

(6.38) (/ci - l) 3 (3a/cj - 1) = 0. 

Note that if we formulate the problem in terms of the extrapolation boundary conditions 
(2.9), then (6.38) follows by inserting (6.33) into (2.9a). 

Example 6.1. As our first example we let a = 0 and (6.38) becomes 

(6.39) (/ci - l) 3 = 0. 

The repeated roots of (6.39) are 

(6.40) /ci = 1 

and from (6.36) the corresponding eigenvalue is A = 0. Recall that the roots of the cubic 
equation (6.36) are denoted by «i,/C2, and /C2. By inserting A = 0 into (6.36), one finds 
that the roots are 

(6.41a, b,c) /ci = 1, k 2 = (5 - \/33)/2 « -0.372, /c 3 = (5 4- V / 33)/2 « 5.372. 
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It is obvious that inequality (6.32) is not satisfied. Therefore, if a = 0 the left boundary 
condition does not introduce any boundary condition dependent eigenvalues. 

Example 6.2. Although the boundary condition (6.34) with (6.35) and a = 0 does not 
introduce any boundary condition dependent eigenvalues, they are introduced for a range 
of nonzero values of the parameter a. The roots of (6.38) are 

(6.42) /ci = 1, 1, 1, l/(3a) (a ± 0). 

Are there any values of a for which A = 0 is a boundary condition dependent eigenvalue? 
For A = 0, the roots of (6.36) are 

(6.43) k = 1, /c = (5 — \/33)/2 » —0.372, k = (5 + >/33)/2 « 5.372. 

It is obvious that inequality (6.32) is not satisfied if k x = 1 which excludes the repeated 
roots in (6.42). Inequality (6.32) is satisfied if 

(6.44a, b,c) = (5 — \/33) / 2 ~ — 0.372, /C2 — 1, /C3 = (5 + \/33)/2 ~ 5.3/2. 

In order to have a boundary condition dependent eigenvalue A = 0 we must choose or in 
(6.42), i.e., 

(6.45) /ci = l/(3a) 

so that is (6.44a). By equating (6.45) and (6.44a), one obtains 

(6.46) ck = — (5 + \/33)/12 ~ — 0.895. 

Consequently, for this value of a the boundary condition (6.34) with (6.35) will introduce 
the boundary condition dependent eigenvalue A = 0. 

If one carries out a GKS normal mode analysis for the semi-discrete approximation 
with (2.16) and the numerical boundary condition (6.34) with (6.35), one finds that the 
semi-discrete approximation is unstable for ct < or where dr is defined by (6.46). 

As explicit formula for A in terms of the parameter a is found by substituting (6.45) 
into (6.37): 

(6.47) A = (3ar - l)(-18a 2 - 15a + l)/(54a 2 ). 

If one chooses a particular value of a , the corresponding value of A is a boundary condition 
dependent eigenvalue if and only if the corresponding /c’s satisfy inequality (6.32). 
Example 6.3. Is it possible to choose or so that there is an eigenvalue 

(6.48) A = -4/3 

which falls on the oval shaped circulant spectrum where the oval crosses the real axis (see 
Fig. 5.5a)? This value of A is an eigenvalue of the circulant matrix (1.6) with (2.16) as one 
can easily verify from the formula (2.17) for 9 = n. (The corresponding /c for the circulant 
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case is k = — 1 .) By substituting A = —4/3 into the cubic equation (6.36) one obtains the 
roots 

(6.49a, b,c) k = - 1 , k = (7 - v^l)/2 « .298, k = (7 + V5I)/2 « 6.70. 

Inequality (6.32) is satisfied if and only if 

(6.50a, b,c) «i = (7 - v / 41)/2 « .298, k 2 = -1, k 3 = (7 + v^l)/2 « 6.70 

If a in (6.45) is chosen so that K\ is (6.50a), one obtains 

(6.51) a = l/(3«i) = (7 + V5T)/12 « 1.117. 

Hence A = —4/3 is a point of the boundary condition dependent spectrum for the value of 
a given by (6.51). 

6.4.2. Right boundary dependent spectrum. In the analysis for the boundary 
condition dependent part of the spectrum, the left and right boundary conditions are 
uncoupled. It is convenient to have a single algorithm for both the left and right boundary 
condition problems. This can be done by switching the matrix so that the right and left 
boundary conditions are interchanged. For example, if the matrix (1.4) is switched as 
described in Section 4.2 one obtains the the matrix: 

'On Cl 2 C 13 C 14 
C 2 1 C 2 2 C 23 C 24 

CL— 2 d — 1 Qo O 

(6.52) A = .... 

O a - 2 a_i do di 

a _2 O—l Oo CL \ 
b\4 bi3 bi2 fen 

where aj = a-j . For this matrix p = 2 and q = 1 and the left boundary difference equations 
are given by 

(6.53a) Cll^l + C\2<f>2 + Cl 3<^3 + C\\(f >4 = A<^i 

(6.53b) C2X<j>l + C 22 <f>2 + C 2 3<^3 + C 2 4^4 = A^ 2 . 

For the matrix (6.52), the algebraic equation (2.6) is 

(6.54) A = d_ 2 /c -2 + d_iK -1 + do + a,\K. 

There are two boundary conditions (6.53) and following the algorithm for the boundary 
condition dependent spectrum we seek a solution of the form (4.5), i.e., 

(6.55) 4>j = /3ik\ +02*2- 
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From (6.53a,b) one obtains 
(6.56a) 

Ci i(/?i K\ + 02^2 ) 4" c i 2 (/^l “b /^2 ^*-2 ) ~b *b 02^2) "b ^14^01 ^1 "b 02^2 ) = '^(/^l K l'b 02 ^2 ) 

(6.56b) 

C2l(/?1«1 + 02*2) + C22(fll + 02^1) + c 23(/^l H~ /?2 ) +C24(/?1 +^2^2) = ^(^1 K 1 + #2 K 2)- 

Since /cj and k 2 must each give the same A, we use (6.54) to obtain two additional equations 
(6.56c) A = d_2«J" 2 + a-i^i -1 + <*o + ai/ci 

(6.56d) A = a_ 2 «^ 2 + d-iK^ 1 + d 0 + Si«2- 

The system of (four) equations (6.56) has (five) unknowns 0\, 02, «i , K 2 , A. Since the equa- 
tions (6.56a,b) are linear and homogeneous in 0\ and 0 2 and (6.56) are algebraic equations, 
they can be reduced to a single algebraic equation in a single variable (we choose k 2 ), i.e., 

( 6 . 57 ) P(k 2 ) = 0 . 

The roots of (6.57) are all possible values of k 2 which may introduce boundary condition 
dependent eigenvalues. Each k 2 is tested as follows: Substitute k — k 2 into (6.o4). Solve 

the new equation for A(k 2 )- Substitute A(k 2 ) back into (6.54) and solve the resulting 

equation for /Cj, «2j K 3- If the K S satisfy (3.1) with (3.4), i.e., 

(6.58) |«i|<M<M 

then A(« 2 ) is a boundary condition dependent eigenvalue, otherwise it is not part of the 
asymptotic spectrum. 

If the bandwidth of the matrix is large and the number of rows in the boundary condi- 
tion matrix B becomes large, the algebraic reduction of the system, e.g., (6.57), becomes 
hopeless even for good symbolic systems such as MACSYMA. In these cases we recommend 
the less reliable but more efficient algorithm of Section 4.3. 

An alternative solution procedure for (6.56) is the following. The system of equations, 
e.g., (6.56) is symmetric in the /c’s so the algebra can be simplified by the introduction of 
the elementary symmetric functions. As an example we consider the relatively simple right 
boundary problem for the matrix (1.7) and introduce the elementary symmetric functions 
[ 11 ] 

(6.59) y = K 1 + AC2 

and 

(6.60) x = kik 2 . 
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It can be shown that the system (6.56) reduces to 


(6.61) 


( \2 5 i 2 cj 3 a_iai 4 . ^13®— i (di ci 2 )ai , 3 

— Ci 3 (- Yx° + X +[ Tj + = \ x 


a - 2 


— 2 


u ._ 2 

+ [C13 - 


a _ 2 


(aj — ci 2 )a_i 2 , \ - n 

— ]x + (a 0 — ciijx — a_ 2 = 0 


a _ 2 


where cn = —3/3, cj 2 = 3/3 — 1/2, and C13 = — /3 and 
(6.62) y = (dix 2 — a_ix)/a_2. 

It is not practical to proceed analytically so for a particular set of matrix coefficients we 
solve (6.61) numerically for x and proceed as in the general algorithm 4.2. 

7. Eigenvector scaling. If conventional numerical eigenvalue packages, e.g., IMSL, 
EISPACK, etc., are used to compute the spectra for non-normal Toeplitz matrices, large 
errors in the spectra may be encountered for matrices of relatively low order. For example 
if we choose the matrix 

(7.1) [-1/6,1, -1/2, -1/3] 


and use the MATLAB routine eig( A) on an IRIS workstation to compute the eigenvalues, 
we obtain the sequence of spectra (for increasing matrix order) shown in Figs. 7.1a through 
7. Id. The Toeplitz matrix defined by (7.1) is the transpose of the Toeplitz matrix defined 
by (2.16). 

If one computes the numerical spectra of the Toeplitz matrix (2.16), e.g., see Fig. 5.1, 
for J = 80 and 160, they would fall on the asymptotic spectra. In general, if one computes 
the numerical spectra of a nonnormal matrix A and its transpose A T (for sufficiently 
large J) they will differ in spite of the fact that the analytical spectra are identical. This 
phenomenon is discussed at the end of this section. 

A non-normal Toeplitz matrix ( e.g., (7.1)) can have the property that the spectrum is 
very sensitive to perturbations of the matrix elements. For such matrices, Trefethen [12] 
and Reichel and Trefethen [8] suggest that an eigenvalue analysis may lead to incorrect 
conclusions and it is more meaningful to analyze pseudo-eigenvalues. The erroneous eigen- 
values plotted in Fig. 7.1 axe a manifestation of the pseudo-spectra of the matrix (7.1) 
for increasing J. The numerical error can be reduced by using higher numerical preci- 
sion. However, for any numerical precision the qualitative behavior illustrated in Fig. 7.1 
will always appear for sufficiently large J. Qualitatively, the numerical Toeplitz spectrum 
approaches the spectrum of its Toeplitz circulant cousin as J — » 00. 

The difficulty in numerically computing the eigenvalues of a non-normal Toeplitz matrix 
is related to the exponential character of the eigenvectors. We illustrate this with an 
example by comparing the eigenvectors of the matrix (2.16) with the eigenvectors of its 
transpose (7.1). Since (2.16) is a Toeplitz matrix, the asymptotic spectrum is independent 
of the boundary conditions and the k roots satisfy (6.23). In this example the moduli 
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FIGURE 7.1. Numerically computed spectra ( o symbol) for the unsealed Toeplitz ma- 
trix [— 1/6, 1, —1/2, —1/3]. Asymptotic Toeplitz spectra and asymptotic circulant spectra 
indicated by “solid” and dashed lines, respectively. 

of the roots of equal modulus (4.1), |/c|, are only weakly dependent on the angle tp( as 
illustrated in Fig. 7.2. From (4.4) with coefficients given by (2.16) and tp = tt/2 one 
obtains |k| = 1/a/ 3. Hence for the roots of equal modulus, one finds 

(7.2) |/ci| = |/«2 1 = |«| « l/\/3 

and consequently the moduli of the eigenvector elements behave as 

(7.3) « (1/V3)> 

i.e., they exponentially decrease with increasing j. 

For the transpose A T , defined by (7.1), one has 

(7.4) |*i| < |*a| = |*3| 

where the check symbol indicates that the k’s in (6.23) and (7.4) are not the same. In fact, 
the /c’s of (6.23) and (7.4) are reciprocals. Consequently, for A T one has 

(7.5) ]«2 1 = |«3 1 = |k| « >/3 
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Figure 7.2. Modulus of k vs i/t for (a) matrix (2.16) unsealed, u, and scaled, s, and (b) 
the transpose matrix (7.1) unsealed, u, and scaled, s. 


and the moduli of the eigenvector elements 

(7.6) l^il ~ (v^y 

grow exponentially with increasing j. 

A simple scaling attenuates the exponential behavior of the eigenvectors. Let A denote 
an arbitrary banded Toeplitz matrix. If the /c’s of equal modulus (3.3) are not on the unit 
circle, then either A or A T will have exponentially growing eigenvectors. Assume that A 
has /c’s of equal modulus (4.1) which are outside the unit circle. Define k to be 

(7.7) k = mean (|/c|). 

0<t/v<7T 

For the matrix A we rescale the eigenvectors by 

m h = If 

The scaling is effectively a normalization of the /c’s so their moduli are approximately equal 
to unity. 

The scaling similarity transformation is 


‘ 1/k 

l/k 2 O 

(7.9) S _1 = 

O l/k 1 - 1 

l/k J 
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FIGURE 7.3. Numerically computed spectra (o symbol) for (a) the unsealed Toeplitz 
matrix [-1/6, 1, —1/2, —1/3] and (b) the scaled matrix. Asymptotic Toeplitz spectra and 
asymptotic circulant spectra indicated by “solid” and dashed lines respectively. 

Hence the rescaled eigenvalue problem is 

(7.10) (S -1 AS)S -1 ^ = XS~ l <f>. 

If a.ij = dj-i are the nonzero elements of the Toeplitz matrix A, the elements a t} of the 
rescaled matrix S -1 AS are 

(7.11) aij = dj-ik j -\ 

To numerically compute the eigenvalues of the rescaled eigenvalue problem, one simply 
replaces the matrix elements by (7.11). One should not make the numerical matrix multi- 
plies indicated by S -1 AS because of potentially large roundoff errors. The implementation 
of algorithm 4.1 in the Appendix computes k(rp) which can be used to compute (7.7). 

In Fig. 7.3 we compare the computed (MATLAB) eigenvalues of the scaled and unsealed 
Toeplitz matrix defined by (7.1). One can predict a bound on the numerical eigenvalue 
distribution of a Toeplitz matrix for large J. In Fig. 7.4a the outer dashed line is the 
circulant cousin spectrum of the Toeplitz matrix defined by (7.1) and the inner dashed line 
is the circulant cousin of the corresponding scaled Toeplitz matrix. The solid line is the 
asymptotic Toeplitz spectrum. For the unsealed matrix the outer dashed line of Fig. 7.4a 
is the bound on the numerically computed eigenvalues for large J. For the scaled matrix, 
the inner dashed line is the asymptotic bound on the numerically computed eigenvalues. 
For the scaled matrix the circulant spectrum is tightly K wound around the asymptotic 
Toeplitz spectrum and one can expect an accurately computed spectrum for large J . These 
bounds are illustrated by the numerical spectra shown in Figs. 7.3. 
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FIGURE 7.4. Asymptotic circulant spectra for unsealed (“outer” dashed lines) and scaled 
(“inner” dashed lines) matrices. Asymptotic Toeplitz spectra are indicated by “solid” lines. 

As a second example, we plot in Fig. 7.4b the spectra of the unsealed and scaled circulant 
cousins of the pentadiagonal Toeplitz matrix defined by 

(7.12) [-1.5, 16, -15,0, .51. 

The “solid line”, i.e., the propeller, is the asymptotic Toeplitz spectrum. Numerically 
computed (MATLAB) spectra of the Toeplitz matrix and its transpose are plotted in Fig. 
7.5a,b for J = 400. The eigenvalues for both A and A T are erroneous for large J . The 
inner dashed curve of Fig. 7.4b is the asymptotic bound on the numerical eigenvalue 
distribution of the scaled matrix for large J. Numerical eigenvalue computations for the 
scaled matrix are shown in Fig. 7.5c,d. 

One can also predict how “well” the scaling will work. If a graph of |/c| vs xpe (see 
algorithm 4.1) for the unsealed matrix is a single “flat” band as in Fig. 7.2 then the 
scaling will yield a tight bound for the circulant cousin around the asymptotic Toeplitz 
spectrum as illustrated in Fig. 7.4a. On the other hand, if |k| vs tp£ forms multiple and/or 
nonuniform bands as illustrated in Fig. 7.6, then one has multiple scales and a simple 
scaling will not yield a closely wound circulant cousin spectrum as illustrated in Fig. 7.4b. 

If one starts with a quasi- Toeplitz matrix, then the appropriate scaling is the same as 
for its pure Toeplitz cousin. 

8. The generalized eigenvalue problem. If Pade type difference approximations 
are used to approximate differential equations the associated stability analysis leads to 
the generalized eigenvalue problem. In addition, the analysis of implicit time integration 
schemes leads to a generalized eigenvalue problem. Let A and B represent J x J banded 
Toeplitz matrices. The generalized eigenvalue problem is defined by 

(8.1) A<t> = \&<f> 
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(c) eig(scaled( A)), J = 400 



( d ) eig(scaled( A T )), J = 1000 


FIGURE 7.5. Numerically computed spectra (o symbols) for (a) A = [—1.5, 16, —15 , 0, .5], 
(b) A t , (c) scaled A, J = 400, and (d) scaled A T , J = 1000. 

where <f> is the eigenvector and A is the eigenvalue. Here we assume that the matrix A is 
defined by the sequence [a,-, ~pa < i < ?a] and B is defined by the sequence [£>,, — p# 
i < qs) where pA,qA,PB,QB are positive integers. The indexing of the eigenvector <f> is 
given by 

(8.2) <}> T = [^i, <f>2, • ■ • , <f>j]- 

The eigenvalue problem (8.1) is equivalent to the set of homogeneous difference equations 


^ ^ <j m (f>j+ m — A ^ ^ b m <f> jA-mi j 1? 2, • ■ ,7 


m=-p A 


m=—pB 


with p = max[p^,ps] homogeneous boundary conditions at the left boundary and q = 
maxfgA) ?b] homogeneous boundary conditions at the right boundary : 


(8.4a, b) 


<t>-m = 0, m = 0, ,p- 1, 


<f>j+m =0, m = 1,2, • ,q. 


* 
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FIGURE 7.6. Modulus of k vs ifit for (a) matrix [—1.5, 16, —15 , 0, .5] unsealed, u, and 
scaled, s, and (b) the transpose matrix unsealed, u, and scaled, s. 

Equation (8.3) is a (p + q)- th order difference equation for (f>j and we look for a solution 
of the form 

(8.5) <f>j = K j 

where k is a complex constant. Substitution of (8.5) into (8.3) yields 

7fl 9A 

(8.6) A 22 b m K m = 22 a m« m - 

m=-pB m=—pA 

This is an algebraic equation of degree p + q in the unknown k and we assume the roots 
are distinct. The general solution of (8.3) is 

p+q 

(8.7) <j)j = 22 ft mK m 

771 — 1 

where the /? m ’s are arbitrary constants. The constants /3 m are determined by inserting 

(8.7) into the boundary conditions (8.4). If the bandwidth is three, i.e., p = q = 1, the 
eigenvalue problem for the pure Toeplitz problem can be solved analytically [15]. However, 
if the bandwidth is greater than three, the problem is analytically intractable. 

The eigenvalue problem for the circulant cousin of (8.1) is (8.3) with periodic boundary 
conditions (2.12). Insertion of (8.5) into (2.12) yields 

( 8 . 8 ) k j = 1 . 
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FIGURE 8.1. Spectra for the generalized Toeplitz eigenvalue problem (8.1), (a) A = 
[2, 0,-1] and B = [-2,2,1] and (b) A = [.38, .13, -.43, .15,^25, -.06, .27, -.02, -.26] and 
B = [.28*, .38*, .50*, .07*. .95*. .99*, .49*, .27*, .09*']. 


Hence the k’s are the J roots of unity given by (2.14) and the circulant spectrum is obtained 
by substituting (2.14) into (8.6): 

gB 

(8.9) A, b m e im9 ‘ = £ a m e im9t . 

m——p b m=—pB 

The partition of the spectrum for the asymptotic Toeplitz spectrum is the same (Section 
3) as for the standard eigenvalue problem. The algorithm follows that of Section 4.1. To 
construct the polynomial for k’s of equal modulus we use (4.1). Since K a and k\, (defined 
by 4.1) must each give the same value of A when inserted in (8.6), we can eliminate A 
and obtain a polynomial in k with coefficients which are functions of the Toeplitz matrix 
coefficients and vt ■ In particular, 

?B ?A 

(8.10a) A ( Ka ) Y b m k m e im ^ = Y 

m = —p b m — ~P A 

qB gA 

(8.10b) A(k 6 ) Y b m k m e~ im ^ = Y 

m=—pB m=—PA 

and since 

(8.11) A(k 0 ) = A (kj) 
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one obtains, after some algebraic simplification, 


qs 

(8.12) ^2 a mb n K ( ' m+n ^ sin[(m - n)ipt] = 0. 

m=—pA n= p b 

The algorithm follows that of (i) and (ii) in Section 4 . 1 . The implementation of the gen- 
eralized eigenvalue algorithm is a straightforward modification of the standard eigenvalue 
algorithm (Appendix). 

As an example of the k polynominal (8.12), if we choose A and B to be the quadridiago- 
nal matrices represented by the sequences [a_2, a~i, uo> Gi] and [6_2, 6— i , &o, &i], respectively, 
the equation for the k ' s of equal modulus is 

( 8 . 13 ) 

(ai&o — aobi)k 4 + 2 (aib-\ — a-ibi) cos(ipe)k 3 + 

[(a 0 b - 1 - a-ibo) + (ai6_ 2 - a_ 2 6 i)( 4 cos 2 (^) - 1 )]k 2 + 

2(a 0 6_2 — a-2bo) cos(ipt)k + (a_ii>_2 — a_2&-i) = 0. 

Note that ( 6 . 4 ) and the “transpose” of ( 6 . 25 ) are easily obtained as special cases of ( 8 . 13 ). 
As examples of the spectra of the generalized eigenvalue problem we choose two examples. 
The spectra for the Toeplitz and the associated circulant Toeplitz matrices are plotted in 
Fig. 8.1 

9. Concluding Remarks. The algorithms presented in Section 4 provide a practi- 
cal accurate method for obtaining the asymptotic spectra of banded Toeplitz and quasi- 
Toeplitz matrices. They are especially useful for non-normal matrices where conventional 
algorithms may give erroneous results. The extension of the algorithms for the generalized 
eigenvalue problem for Toeplitz matrices was presented in Section 8. The asymptotic anal- 
ysis also provides a similarity transformation, Section 7 , which can increase the spectrum 
accuracy if conventional eigenvalue algorithms are employed. 

Acknowledgements. We are grateful to Dennis Jespersen, Robert Schreiber, and Nick 
Trefethen for comments and suggestions which improved our manuscript. 

Appendix. MATLAB Toeplitz eigenvalue algorithm. 

function[teig, spsi, sakap] = toeasy(c,r,nn) m , 

%Find the asymptotic (size approaches infinity) spectrum 
%for a banded Toeplitz matrix with first column V and first row 
%’r’ (only banded part of column and row) nn is a measure of the 
% number of points computed on the spectrum 
diag = c(l); 

Ir = length(r ); 
nr = Ir — 1; 

1 c = length(c)\ 
nc — Ic— 1; 
ct = c( 2 : /c); 
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rt = r( 2 : /r); 
mu — nr : —1 : — nc; 
pcO = [ fliplr(rt ) 0 ct]; 
teig = [1 : (nc + nr) * nn); 
spsi = [1 : (nc + nr) * nn]; 
sakap = [1 : (nc + nr) * nn]; 
eigct = 0; 
epst = 1000 * eps; 
for j = 1 : nn; 
psi = j * pi/(nn + 1); 
mvs = sin(mv * psi); 
pe = pcO. * mvs; 
nroot = nc + nr; 
i/ (pe(2) ~= 0); 

if(abs(pe(l)/pe(2 )) < 1000 * eps); 
pe = pe( 2 : nroot + 1); 
nroot = nroot — 1; 
end; 
end; 

ke = roots(pe); 
for k = 1 : nroot; 
if ke{k) ~= 0; 
kap = ke(k) * exp(i * psi ); 

/am = —(kap. A mv) * pcO.'; 

pc = [//ip/r(rf) lam ct]; 

kc = roots(pc ); 

akap = abs(kap); 

akc = abs(kc); 

skap = sort(akc); 

test = abs((skap(nc) — skap(nc -f 1 ))/(skap(nc) + skap(nc + 1))); 
if test < epstSzskap(nc) > akap * (1 — epst)&cskap(nc) < akap * (1 + epst); 
eigct = eigct + 1; 
teig(eigct) = —(lam — diag ); 
spsi(eigct) = psi; 
sakap(eigct) = akap; 
end; 

end; 

end; 

end; 

teig = te/<7(l : eigct); 
spsi — spsi(l : eigct); 
sakap = sakap( 1 : eigct); 
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